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(f) , Abstract. We study a simple model for the solar dynamo in the framework of the Parker migratory dynamo, 

with a nonlinear dynamo saturation mechanism based on magnetic helicity conservation arguments. We find a 
parameter range in which the model demonstrates a cyclic behaviour with properties similar to that of Parker 
dynamo with the simplest form of algebraic a-quenching. We compare the nonlinear current helicity evolution in 
this model with data for the current helicity evolution obtained during 10 years of observations at the Huairou 
Solar Station of China. On one hand, our simulated data demonstrate behaviour comparable with the observed 
phenomenology, provided that a suitable set of governing dynamo parameters is chosen. On the other hand, the 
l' observational data are shown to be rich enough to reject some other sets of governing parameters. We conclude 

1~ \ that, in spite of the very preliminary state of the observations and the crude nature of the model, the idea of 

using observational data to constrain our ideas concerning magnetic field generation in the framework of the solar 
' dynamo appears promising. 
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I 1. Introduction a (quasi)stationary wave which can be compared with the 

■ - i observed activity cycle, instead of a dynamo wave with an 

The solar activity cycle is widely believed to be connected ,. n . ,., , T . , , 

. J J .... exponentially growing amplitude, m principle, the phe- 

with dynamo action which occurs somewhere inside the so- i r , i i i u j j • 

J . nomenology ot the solar cycle can be reproduced using a 

lar convective zone or even in the overshoot layer. Starting . ... , i r , , 

J ° very primitive a-quenchmg model ot dynamo saturation, 

from the seminal paper of Parker (1955), various dynamo t i i « u 

11 ; ; i with the energy of the dynamo generated magnetic held 

models have been suggested for the solar cycle (see e.g. -,■ ■ . . . . . , , , . . 

aa acnievmg approximate cquipartition with tnc kinetic cn- 

Riidiger & Brandenburg 1995: Soha et al. 1998; Tavakol et t , , • 

° , o ergy ot tne random motions, 
al. 2002; Brooke et al. 2002; also Blackman & Brandenburg 

2003, whose dynamo model of solar cycle also includes A dee P er treatment of solar dynamo saturation re- 
magnetic helicity balance). These models exploit particu- q uires however some ideas concerning the physical pro- 
lar parameterizations for sources of the dynamo activity, cesses that S ive rise to quenching of the generation mech- 
i.e. the a-effect, which in turn is connected with the mean anism - A scenario of dynamo saturation which is now 
hydrodynamic helicity X v of the convective motions, and widel y discussed is connected with the concept of magnetic 
acts in conjunction with the nonuniform rotation il. If the helicity. The point is that the weakest link in the dynamo 
dynamo action is strong enough, a dynamo wave propa- self-excitation chain, i.e. a, is a pseudoscalar quantity and 
gating somewhere inside the convective shell is excited. It cannot be directly connected with the magnetic energy, 
is necessary to include some saturation mechanism to get which is a scalar (not pseudoscalar) quantity. A magnetic 

helicity x m can however be introduced to describe the level 

Send offprint requests to: I. Rogachevskii of magnetic field mirror- asymmetry and this quantity can 
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be associated with the magnetic part of a, i.e. a m , which 
is thought to be responsible for a-quenching. 

The magnetic helicity x m is an integral of motion for 
the ideal MHD equation, similar to the hydrodynamic he- 
licity x v which is conserved in the hydrodynamical case. 
During the solar activity cycle, magnetic helicity is redis- 
tributed between the large and small scale magnetic field. 
Based on this concept, a governing equation for a m has 
been suggested (Kleeorin & Ruzmaikin 1982; Kleeorin & 
Rogachevskii 1999). Together with the mean-field dynamo 
equations, this equation has solutions in form of a prop- 
agating steady dynamo wave (Kleeorin et al. 1994, 1995; 
Covas et al. 1998; Blackman & Brandenburg 2002). 

For a long time, it was impossible to observe either the 
magnetic helicity x m or the hydrodynamic helicity \ v an d 
these values, crucial for dynamo theory, were taken from 
theoretical estimates only. In last decade, basic progress 
here has been achieved and the first observations of mag- 
netic helicity in active regions on the solar surface have 
been obtained (Pevtsov et al. 1994, 1995; Zhang & Bao 
1998, 1999; Canfield & Pevtsov 1998; Longcope et al. 
1998). It is possible to some extent to isolate a latitudinal 
distribution of magnetic helicity averaged over a solar cy- 
cle (Zhang et al. 2002) as well as to follow the temporal 
evolution of magnetic helicity averaged over latitude. The 
obvious aim now is to confront predictions of dynamo the- 
ory concerning the latitudinal distribution of magnetic he- 
licity and its evolution during a solar cycle with the corre- 
sponding observational data; this is the aim of the present 
paper. When carrying out our investigation, we take into 
account that the available data concerning magnetic helic- 
ity of the solar magnetic field are still quite uncertain, and 
it would be unrealistic to expect that more or less fine de- 
tails can be isolated using this data. Correspondingly, we 
restrict ourself to a very crude theoretical model, that we 
confront with observations. Specifically, we simplify the 
mean-field dynamo equations at the level of the Parker 
migratory dynamo equations, and include the algebraic a- 
quenching and the dynamic a-quenching associated with 
magnetic helicity evolution as the only saturation mech- 
anisms. Both hypotheses are obvious simplifications and 
there is no problem in principle in including many more 
realistic features into our dynamo model. However we con- 
sider that to be a topic for further work. 

The other point to be clarified from the very beginning 
is the following. Magnetic helicity can be understood as a 
measure of the linkage of magnetic lines and it is necessary 
to reconstruct the complete 3D magnetic field structure 
to deduce this helicity from observations. Clearly, this is 
a very complicated observational problem and various in- 
termediate quantities such as current helicity (b- (Vxb)), 
i.e. the linkage between electric current lines, are used to 
this end (where b represents magnetic fluctuations). These 
quantities are useful in theoretical studies of dynamo sat- 
uration also and we use them below. Our work needs a 
clear distinction between such concepts as a-effect and 
the corresponding helicity, and between hclicities of total 
magnetic field, large-scale magnetic field and small-scale 



magnetic field; these distinctions can be neglected to some 
extent in other areas of dynamo theory. 

2. Magnetic and current helicity data obtained at 
the Huairou Solar Observing Station 

The averaged value of the small-scale magnetic helicity, 
i.e. (a • b), evidently would be a convenient quantity to 
confront with a dynamo saturation scenario based on a 
magnetic helicity conservation argument. In practice how- 
ever, the vector potential a, being a non-gauge invariant 
quantity is inconvenient observationally, and it is the cur- 
rent helicity (6 • (V x b)) = (b x (V x b) x ) + (6j,(V x b) y ) + 
(i>z(V x b) z ) which can be extracted from the observations 
(here x, y, z are local cartesian coordinates connected with 
a point on the solar surface and the z-axis is normal to the 
surface). The observations are restricted to active regions 
on the solar surface and we obtain information concerning 
the surface magnetic field and helicity only. Monitoring of 
solar active regions while they are passing near to the cen- 
tral meridian of the solar disc enables observers to deter- 
mine the full magnetic field vector. The observed magnetic 
field is subjected to further analysis to obtain the value 
V x 6. Because it is calculated from the surface magnetic 
field distribution, the only electric current component that 
can be calculated is (V x b) 2 . As a consequence of these 
restrictions, the observable quantity is 

H c = (MV x b) z ) . (1) 

Because the surface magnetic field is almost force- 
free, it is also useful to consider the magnetic field twist 
as (Woltjer, 1958) which is defined as the proportional- 
ity coefficient between magnetic field and electric current 
(j = Offb). The observational restrictions discussed above 
imply that the observational equivalent of twist is 

as=jz/b z . (2) 

The observational data used in our analysis were ob- 
tained at the Huairou Solar Observing station of the 
National Astronomical Observatories of China. The mag- 
netographic instrument based on the Fel 5324 A spectral 
line determines the magnetic field values at the photo- 
spheric level. The data are obtained from a CCD camera 
with 512 x 512 pixels over the whole magnetogram, whose 
entire size is comparable with the size of an active re- 
gion, as well as with the depth of the solar convective zone 
(about 2 x 10 s m). However, because of the observational 
technique, the line-of-sight field component b z can be de- 
termined with a much higher precision than the transverse 
components (b x and by). There are a number of other ob- 
servational difficulties such as in resolving the so-called 
'180° ambiguity' in the direction of the transverse field 
etc. The observational technique is described in detail by 
Wang et al. (1996), see also Abramenko et al. (1996). 

An observational programme to reveal the values of the 
twist and the current helicity density over the solar surface 
requires a systematic approach, both to the monitoring of 
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Table 1. The first column gives the central latitude 9 = 
90° — 9 of the data bin, with the averaging interval in 
brackets, the twist (ag) is measured in units of 10 _8 m _1 , 
the current helicity (H c ) in units of 10~ 3 G 2 m -1 , and N is 
the number of active regions involved in the analysis. The 
errors correspond to the 95% confidence level, i.e. about 
two standard deviations. 



e 


(ag) 


(He) 


N 


28 (24 - 32) 


-0.4 ± 1.2 


-1.6 ± 1.7 


18 


20 (16 - 24) 


-0.9 ±0.8 


-0.9 ±0.4 


51 


14 (12 - 16) 


-1.7 ± 1.3 


-0.6 ±0.4 


34 


10 (8 - 12) 


-2.2 ±0.6 


-0.4 ±0.2 


49 


4 (0 - 8) 


-1.9 ±0.8 


-0.6 ±0.2 


44 


-4 (-8 -0) 


0.3 ±0.7 


0.7 ±0.5 


31 


-10 (-12 - -8) 


1.2 ±0.7 


0.7 ±0.4 


59 


-14 (-16 - -12) 


0.9 ±0.7 


0.9 ±0.7 


46 


-20 (-24 - -16) 


1.0 ±0.8 


0.4 ±0.2 


68 


-28 (-32 - -24) 


1.6 ± 1.7 


0.5 ±0.9 
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Table 2. The data of Table 1 binned by hemisphere and 
year of observation. 



T 




m 


N 


North 


1988-89 


-1.1 ±0.8 


-1.0 ±0.5 


50 


1990-91 


-1.0 ±0.7 


-1.0 ±0.5 


61 


1992-93 


-2.1 ±0.7 


-0.7 ±0.3 


45 


1994-95 


-2.6 ±0.9 


-0.3 ±0.1 


34 


1996-97 


-1.2 ± 1.0 


-0.2 ±0.2 


9 


South 


1988-89 


1.0 ± 1.2 


0.2 ±0.3 


38 


1990-91 


0.9 ±0.7 


0.8 ±0.6 


65 


1992-93 


1.2 ±0.5 


0.9 ±0.3 


77 


1994-95 


0.7 ±0.9 


0.1 ±0.1 


35 


1996-97 


0.3 ±2.0 


0.2 ±0.3 


8 



magnetic fields in active regions and to the data reduc- 
tion, in order to reduce the impact of noise. This work 
has been carried out by a number of research groups (e.g. 
Seehafer 1990; Pevtsov et al. 1994; Rust & Kumar 1996; 
Abramenko et al. 1997; Bao & Zhang 1998; Kuzanyan et 
al. 2000). While this work is still in progress, the largest 
systematic data-set of active regions presently available 
consists of 422 active regions over the 10 years 1988-1996 
(Bao & Zhang 1998). We use averages and confidence in- 
tervals calculated from these data (Table 1). Note that 
both the averaged quantities of twist and current helicity 
density are positive/negative over Southern/Northern so- 
lar hemispheres respectively, and thus obey the so-called 
hemispheric rule (Table 2). 

Observations at Huairou Observing Station also give 
dopplergrams of velocity fields v over active regions from 
the Fel 5324 A spectral line in the photosphere and from 
the Hp line in the chromosphere. The processing of these 
data could provide values of v z b z (which is related to the 
so-called cross-helicity) in the foreseeable future. However, 



at present this quantity is not available for statistical stud- 
ies. 

3. The dynamo model 

We describe the solar dynamo by means of the mean-field 
equation (e.g. Moffatt 1978; Parker 1979; Krause & Radler 
1980; Zeldovich et al. 1983) which in general form is 



dB 

~dt 



Vx(VxB + £ — 77o VxB) 



(3) 



where V is a mean velocity (e.g. the differential rotation), 
770 is the magnetic diffusion due to the electrical conduc- 
tivity of the fluid, £ = (u x b) is the mean electromotive 
force, u and b are fluctuations of the velocity and magnetic 
field respectively and angular brackets denote averaging 
over an ensemble of fluctuations. The electromotive force 
£ can be separated into several contributions, which in- 
clude the a-effect, turbulent magnetic diffusivity r\ and 
other terms such as the magnetic turbulent diamagnetic 
effect. For now, we restrict ourself to the two first terms, 
and consider a and 77 to be isotropic quantities. We take 
the turbulent diffusivity as a prescribed quantity and take 
into account the nonlinearity of the a-effect only, i.e. we 
use the parameterization 



£ = a(B)B - r/VxB 



(4) 



where a depends in principle on the entire evolution of the 
magnetic field, rather on its value in a given instant. This 
dependence is described by an evolution equation (such as 
Eqs. © and JT3JI below.) 

Using spherical coordinates r,9,<fi, we represent an ax- 
isymmetric mean magnetic field as B = B c / 1 e^ + 'Vx(Ae ( j > ). 
Following Parker (1955) we consider dynamo action in a 
thin convective shell, average A and B^ over the depth of 
the convective shell and consider these quantities as func- 
tions of colatitude 9 only. Then we neglect the curvature 
of the convective shell and replace it by a flat slab to get 
the following equations (we drop the suffix on for the 
sake of brevity) 



dB 

~di 
OA 
~dt 



B A d 2 B 

a 2 a 2 „ 

aB + —r- - pi A , 



fi 2 B 



(5) 
(6) 



(see Appendix A) . Here we measure lengths in units of the 
solar radius and time in units of a diffusion time based on 
the solar radius and turbulent magnetic diffusivity. The 
terms —^i 2 B and —/j, 2 A represent the role of turbulent dif- 
fusive losses in the radial direction - the value /1 = 3 corre- 
sponds to a convective zone with a thickness of about 1 /3 
of the solar radius, g — dil/dr is the radial shear of differ- 
ential rotation. We neglect any latitudinal dependence of 
the rotation curve as well as the link between poloidal and 
toroidal magnetic field via the a-effect, and so consider a 
simple aw-dynamo. a and g are normalized with respect 
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to their maximal values and incorporated into the dimen- 
sionless dynamo number D, which gives the intensity of 
the dynamo action (see Sokoloff et al. 1995 for mathemat- 
ical details of the derivation of Eqs. © from Eqs. J3J), 

61). 

These equations are obviously oversimplified. Starting 
from the fundamental paper of Parker (1955) they can 
be used however to reproduce basic qualitative features of 
solar and stellar activity and appear to be viable for this 
purpose. Taking into account the nature of the approach, 
we use the simplest profiles of dynamo generators com- 
patible with symmetry requirements, i.e. a{9) — cos 9 and 
<? = !• 

The points 9 = and 9 = 180° correspond to North 
and South poles respectively. We take here zero boundary 
conditions for A and B. (Because we neglect the convective 
shell curvature, these boundary conditions necessarily are 
approximate.) In principle, we could use a slightly more 
elaborate version of Eqs. I® which take into account 
some curvature effects, so the diffusion term becomes for- 
mally singular at the poles and a more realistic finiteness 
condition can be exploited (see Galitski & Sokoloff 1999). 
Magnetic helicity data are available for middle latitudes 
and the equatorial region (—30° < 9 < 30° where the lat- 
itude 6 = 90° — 9 and 9 = corresponds to the equator) 
only and so here we are not very interested in details of dy- 
namo wave behaviour near to the poles. We keep a factor 
sin 9 in Eq. J5J which reflects the fact that the length of 
the parallels 9 = const vanishes at the poles (Kuzanyan & 
Sokoloff 1995). Neglecting this term results in an unphys- 
ical coupling between the dynamo wave behaviour near to 
the pole and near to the equator. 

We are interested in dynamo waves propagating from 
middle solar latitudes towards the equator. This corre- 
sponds to negative dynamo numbers provided a is chosen 
to be positive in the Northern hemisphere and g is positive 
near to the solar equator. According to various models, the 
ranges of \D\ « 10 3 — 10 6 can be considered as realistic 
for the solar case. 

We nondimensionalize the dynamo equations by mea- 
suring length in units of the solar radius i?, time in units 
of the turbulent magnetic diffusion time R 2 /ij, and the 
differential rotation Sfl in units of the maximal value of 
fi. a is measured in units of the maximum value of the 
hydrodynamic part of the a-effect. 

It is convenient to present the dynamo number as 
D = R Rlu, where R a = aRj-q ~ 1 — 200 and R w = 
SflR 2 jr\ ~ (1—4) x 10 3 (where typical values of parameters 
have been used for these estimates) represent the contribu- 
tions of the a-effect and differential rotation, respectively. 
We use the equipartition magnetic field = u^/Aitp as 
the unit of magnetic field. The vector potential of the 
poloidal field A is measured in units of R a RB*, the density 
p normalized to its value at the bottom of the convective 
zone, and the basic scale of the turbulent motions I and 
turbulent velocity u at the scale I are measured in units 
of their maximum values through the convective region. 
The magnetic Reynolds number Rm = Iu/t]q is defined 



using these maximal values, rj = lu/3 is an estimate for 
the turbulent diffusivity. 

We stress that all physical ingredients of the model 
vary more-or-less strongly with the depth h below the 
solar surface and we have to use some average quanti- 
ties in the Parker dynamo equations. We use mainly esti- 
mates of governing parameters taken from models of the 
solar convective zone (see, e.g., Spruit 1974 and Baker 
& Temesvary 1966; more modern treatments make lit- 
tle difference to these estimates). In particular, at depth 
h ~ 2 X 10 10 cm, Rm - 2 • 10 9 , u ~ 2 x 10 3 cm s -1 , 
I ~ 8 x 10 9 cm, p ~ 2 x lO" 1 g citT 3 , 77 - 5.3 x 10 12 
cm 2 s _1 . The density stratification scale is estimated here 
as A p = p/\V p\ ~ 6.5 x 10 9 cm and the equipartition 
mean magnetic field B* = 3000 G. In the upper part of 
the convective zone, say at depth h ~ 2 x 10 7 cm, these pa- 
rameters are Rm ~ 10 5 , u ~ 9.4x 10 4 cm s _1 , 1 ~ 2.6x 10 7 
cm, p ~ 4.5 x 10~ 7 g cm -3 , rj ~ 0.8 x 10 12 cm 2 s _1 and 
A p ~ 3.6 x 10 7 cm; the equipartition mean magnetic field is 
£?* = 220 G here. This estimate for the equipartition mag- 
netic field at the base of the convection zone (i?» = 3000 
G) is roughly consistent with the magnetic field strength 
in sunspots (about 1 kG). However obviously it should be 
distinguished from the mean magnetic field at the solar 
surface; a deeper discussion of this distinction is outside 
of the scope of the paper. For the Parker migratory dy- 
namo, the toroidal magnetic field usually dominates and 
below we ignore the poloidal magnetic field when calcu- 
lating the magnetic energy. 

4. The nonlinearities 

A key idea of the dynamo saturation scenario exploited 
below is a splitting of the total a effect into the hydrody- 
namic (a v ) and magnetic (a m ) parts 

a(r, 9)=a v + a m , (7) 

as first suggested by Frisch et al. (1975). We need to pa- 
rameterize both contributions, a v and a m , in terms of the 
magnetic field components and helicities. Two types of ef- 
fect should be taken into account. First of all, the link be- 
tween a-effect and the relevant helicities can be modified 
by the dynamo-generated magnetic field. Correspondingly, 
we introduce quenching functions <f> v (a v = x V( Pv with 
X v = -( T /3)(w • (Vxit))) and 4> m (a m = x c <l> m , and \ c is 
defined below) to obtain a = x v 4>v J rX c< t > rm where r is the 
correlation time of the turbulent velocity field. The second 
problem to be addressed is that magnetic helicity is not 
a very convenient quantity because it involves a gauge- 
noninvariant quantity, i.e. the vector potential. We con- 
nect magnetic helicity with the current helicity (b- (V x b)) , 
by using the approximation of locally homogeneous turbu- 
lent convection (see Kleeorin & Rogachevskii 1999). Then 
we need to obtain a quantity of suitable dimension, and in- 
troduce the density p to obtain the correctly dimensioned 
X c = (r/127rp)(6 • (Vxb)). Thus, 

a(r,e)= X v <t>v+X c <t>m- (8) 
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The quenching functions 
given by 



and (f> m in Eq. JSJl arc 



MB) = (l/7)m m (B)+3L(B)} , 
<t>m(B) = ^[l-arctan(V8B)/V8B] 



(9) 
(10) 



(see Rogachevskii & Kleeorin 2000, 2001), where L(B) = 
1 - 16B 2 + 128S 4 ln(l + l/(85 2 )). Thus <j> v = 1/(4B 2 ) 
and (j> m = 3/(8B 2 ) for S > 1/3; and <f> v = 1 - (48/5)B 2 
and (j) m = 1 - (24/5)S 2 for B < 1/3. Here X v and x c are 
measured in units of the maximal value of the a-effect. 

The function <f> v describes conventional quenching of 
the a effect. A simple form of such a quenching, <p v = 
1/(1 + B 2 ) 1 was introduced long ago (see, e.g. Iroshnikov 
1970). This form is quite close to the more sophisticated 
form presented in Eq. (J5J. The magnetic part a m in- 
cludes two types of nonlinearity: the algebraic quench- 
ing described by the function <\> m (see e.g. Field et al. 
1999; Rogachevskii & Kleeorin 2000, 2001) and the dy- 
namic nonlinearity which is determined by Eq. (|13|) . 

The quenching of the a-effect is caused by the direct 
and indirect modification of the electromotive force by the 
mean magnetic field. The indirect modification of the elec- 
tromotive force is caused by the effect of the mean mag- 
netic field on the velocity fluctuations and on the magnetic 
fluctuations, while the direct modification is due to the ef- 
fect of the mean magnetic field on the cross-helicity (see, 
e.g., Rogachevskii h Kleeorin 2000, 2001). 

We can calculate also the cross-helicity (it • b) which 
may in the future be compared with observational data 



(u • b) = (7 7 /2)[3A- 1 B r + <f> ch {B) (B ■ V)B 2 } 



(11) 



where A" 1 = |V(it 2 )|/(u 2 



fell 



(B) = (2/35B 2 )[(15 



224B 2 )(t> m (2B) + 6L(2B) - 21], and c/> ch (B) = -128/5 for 
B < 1/3, and cf> ch (B) = -3tt/(20V2 B 3 ) for B > 1/3. 
When deriving Eq. {TTJ we used Eqs. (A14)-(A17) and 
(A21) of Rogachevskii & Kleeorin (2001). 

Now we need to average Eq. JHJ over the depth of the 
convective zone. The first term in the averaged equation 
seems to be determined by the values taken at some sort of 
mean position in the convective zone, while the situation 
concerning the second term is much less clear, because 
the density used to calculate x° decreases strongly with 
radius. The clarification of this problem is obviously be- 
yond the Parker approximation; however to address this 
problem as far as it possible here we introduce a phe- 
nomenological parameter a by 



a(6) = x v <f>v + ox c 4>r, 



(12) 



where the helicities and quenching functions are associated 
with some sort of mean position in the convective zone. We 
emphasize that below we consider a as a free parameter 
in the context of the averaging process used to derive the 
Parker equations; probably we can only safely assert that 



For the sake of brevity of notation, we keep in Eq. I|12f> 
the same notation as in Eqs. J7J) and JSJ). We stress that 
now we consider a parameterization for helicities which 
depends on colatitude 6 only. Of course, this is not more 
than a phenomenological description, to be improved in 
more detailed models of the nonlinear solar dynamo. 

The function x°(B) is determined by a dynamical 
equation which follows from the conservation law for mag- 
netic helicity (see Kleeorin & Rogachevskii 1999). A gen- 
eral dimensional form of this equation reads 



dx c 
dt 



r 

T 



9tt t] p. 



{S B + V • *) + kA x c 



(13) 



Here $ = Cx"4'vB 2 l 2 e r /A p is a nonadvective flux of 
the magnetic helicity (here e r is the unit vector in radial 
direction), -kVx c is the diffusive flux of the magnetic he- 
licity (see Kleeorin & Rogachevskii 1999; Kleeorin et al. 
2000; 2002; 2003), and T = l 2 /rj Q is the relaxation time 
of magnetic helicity. Magnetic helicity transport through 
the boundary of a dynamo region is reported by Chae 
(2001) to be observable at the solar surface. We also take 
into account that for an axisymmetric problem the term 
which determines the advective flux of the magnetic helic- 
ity, V • ( Vx c ), vanishes ( V = fi r sin 9 is the differential 
rotation). The parameter C is a numerical coefficient and 
k is of order r\. In principle, these parameters can be cal- 
culated given some model of convection, however we here 
take into account our real level of knowledge and keep 
them as free parameters. Note that in estimating the he- 
licity flux $ we have to include density gradients in the 
radial direction which arc neglected in other parts of the 
analysis. 

The physical meaning of Eq. I|13fl is that the total mag- 
netic helicity is a conserved quantity and if the large-scale 
magnetic helicity grows with magnetic field, the evolu- 
tion of the small-scale helicity should somehow compen- 
sate this growth. Compensation mechanisms include dis- 
sipation and various kinds of transport. 

The dynamical equation (|13(l for the function x c (B) 
in nondimensional form in the context of the Parker mi- 
gratory dynamo reads 



dx'^-i^ 2 , c (2R\ 2 fdAdB d 2 A 
— + (T =( T j 



-aB 2 + 2p 2 AB + CB 2 $ vX v {Q) +k 



d 2 x c 
~d6 2 ~ 



(14) 



(see Appendix A), where we have averaged Eq. (|13f) over 
the depth of the convective zone, so that the averaged 
value of T^ 1 is 



T -l = H -i 



T~ x (r)dr 



Aj R 2 yo 
Hl 2 V 



0.2-0.5 



(15) 



A; is the characteristic scale of the variations /, T(r) = 
(j]/R 2 )(l 2 /r/o) is the nondimensional relaxation time of 
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Table 3. Parameters of activity cycles: B max is the value 
of the dimensionless amplitude of the toroidal magnetic 
field in units of the equipartition field (estimated as 3000 
G at the bottom of the convective zone), T c is the di- 
mensionless cycle length, SP and DP denote singly and 
doubly periodic solutions respectively; 'runs away' means 
that no stable finite amplitude solution was found, al- 
though in some cases there is a long pseudo-stable initial 
phase. The other governing parameters are a = 1, T = 3, 
(2R/1) 2 = 300. For the DP solutions, the amplitude of the 
stronger cycles and the shorter period are given. 



c 


K 


D 


^max 


T c 




0.1 


1 




-10 3 


0.12 


0.37 


SP 


0.1 


1 




-10 4 


5.28 


0.25 


DP 


0.1 


1 


-2 


x 10 4 






runs away 


0.1 


0.1 




-10 3 


0.09 


0.345 


SP 


0.1 


0.1 




-10 4 






runs away 


0.1 


3 




-10 3 


0.13 


0.355 


SP 


0.1 


3 


-3 


x 10 3 


1.47 


0.46 


SP 


0.1 


3 




-10 4 






runs away 


0.01 


1 




-10 3 


0.11 


0.35 


SP 


0.01 


1 


-3 


x 10 3 


1.35 


0.37 


SP 


0.01 


1 




-10 4 


4.50 


0.30 


DP 


0.01 


1 


-2 


x 10 4 






runs away 





1 




-10 4 


4.25 


0.29 


weakly DP 





1 


-2 


x 10 4 






runs away 



magnetic helicity, and the quantities A;, rjo, I in Eq. (|15|l 
are associated with the upper part of the convective zone. 
The parameters C and k can be chosen as C ~ (0.1 — 1) x 
10 _1 and k = 0.1 — 1. The factor 10~ 4 in the coefficient C 
arises from the parameter (l 2 /A p )/R ~ 0.1 (see Appendix 
A). 

5. Results 

We simulated numerically the model of the nonlinear so- 
lar dynamo based on the Parker approximation and con- 
servation of magnetic helicity arguments, as presented in 
previous sections. We found that the model gives a sta- 
ble nonlinear wave-type solution similar to the solar cycle 
phenomenology from the generation threshold of the non- 
linear system (-D n oni,crit ~ —940 with a = k = 1, C = 0.1) 
up to D = — 10 4 . It is quite interesting that a slightly 
stronger generation is needed to get stable nonlinear os- 
cillations with a nonvanishing amplitude than to excite 
the linear dynamo, where Z? crit 1=3 —910. The temporal 
behaviour of nonlinear dynamo waves is quite similar to 
that with the simple algebraic a-quenching, and we give 
here, instead of quite standard plots, details of the cycle 
parameters in Table 3. A feature of the activity cycles illus- 
trated by this Table is a transition from singly to doubly 
periodic solutions with growth of dynamo intensity. 

In the present paper we have concentrated on the dy- 
namics of the current helicity (which is proportional to the 
function x c ) an d its comparison with the observations. At 
the present state of observations we can compare latitudi- 




Fig. 1. The latitude dependence of the time-averaged 
function (x c ) for various parameter sets (see description in 
the text). The observed values of the time-averaged (H c ) 
are shown by filled squares, the error-bars are shown by 
vertical lines. A fitting factor of 23.5, (i.e., (x c )e=28° — 
23. 5 (if c )) has been used - see text. 



nal distributions of \ c averaged over the activity cycle or 
the temporal behaviour of \ c averaged over a hemisphere. 
Such comparisons are presented in Figs. 1-3 for D — — 10 3 , 
a = 1, T = 3, {2R/1) 2 = 300. In Figs. 1 and 3 the param- 
eter k = 0.1. We restrict ourselves to discussing dynamo 
models with singly periodic behaviour (although there are 
hints of a double periodicity in the sunspot record). Two 
types of behaviour are demonstrated. Provided that C is 
negative (magnetic helicity inflow, see Fig. 2b) or small 
and positive (C < 0.1, moderate magnetic helicity out- 
flow, see Fig. 2a) the value of x c is always negative in the 
Northern hemisphere, in accordance with naive theoreti- 
cal expectations as well as the available observations. If C 
is large and positive (strong magnetic helicity outflow), we 
obtain a cycle during which \ c changes sign. Both types of 
behaviour are shown in Fig. 1, for C = 0.01 (solid curve) 
and C — 0.1 (broken). Comparing observational and theo- 
retical results we have to fit the numerical data for \ c with 
its observational equivalent if c . We normalize the data to 
make them equal at 9 = 28° ((x c )e=28° = 23.5(ff c }). 
Because of this fitting procedure, the agreement between 
the observational and numerical data is much better in 
the Northern hemisphere than in the Southern. 

We illustrate the importance of the fitting procedure 
for various parameter sets in Fig. 2 where we give the 
results for C — 0.01 in panel (a) and for C = —0.1 in 
panel (b). The results for k = 0.1 are shown by broken 
curves and those for n = 1 by solid. Both cases represent 
a model which does not exhibit helicity reversals. We ap- 
preciate that the agreement between numerical and obser- 
vational data in the Northern hemisphere can be partially 
attributed to the fitting procedure; however the agreement 
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Fig. 2. The latitude dependence of the time- averaged 
function (x c ) for various parameter sets (see description in 
the text). The observed values of the time-averaged (H c ) 
are shown by filled squares, the error-bars are shown by 
vertical lines. We have used a fitting factor in (a) of 23.5, 
(i.e., (x c )e=28° = 23.5(i? c )) and in panel (b) it is 32. 



obtained looks quite impressive for the primitive models 
considered. A disagreement between the model and ob- 
servations at 6 = 4° can be explained, for example, as a 
result of a non-perfect North-South symmetry in the ob- 
served cycle. We conclude from Fig. 2 that the minimal 
value of (x c ) decreases with k; experience from numerical 
simulations show that it decreases with a as well. Note 
however that the contribution of X c to the magnetic part 
of the a-effect is determined by <r(x c ) and this value is 
more or less stable. 

Moving to the comparison of the temporal helicity 
evolution presented in Fig. 3, we note that the fitting 
procedures now have to be more complicated. The first 
two measured points (1988-89) and (1990-91) with ap- 
proximately equal values of the latitude-averaged (H c ) 
are chosen to be located symmetrically about the mini- 
mum of the function (x c ) for 6 > 0. The distance be- 




3.45 3.5 3.55 3.6 3.65 



Fig. 3. The time dependence of the latitude-averaged 
function ( X c ) for D = -10 3 , a = 1, C = 0.01 and 
k = 0.1. The observed values of the latitude-averaged 
function (H c ) are shown by filled squares (for > - 
Northern hemisphere, lower panel) and filled circles (for 
O < - Southern hemisphere, upper panel), the error-bars 
for (H c ) are shown by vertical lines. The dashed-dottcd 
line indicates the time dependence of the latitude-averaged 
function —6(B 2 ). The filled diamonds in the lower panel 
give the scaled averaged group sunspot numbers, R g - see 
text for further details. 



tween the observational points is (2/ll)(T c /2), where T c 
is the period of oscillations of the toroidal magnetic field. 
The fitting factor is 44.3, determined by the condition 
(x c )t=3.527t d = 44.3(# c ) t= i988, where t d is the turbulent 
diffusion time at the bottom of the convective zone. 

The point is that we can base the fitting procedure on 
a comparison between the simulated value of the toroidal 
magnetic field and some tracer of cyclic solar activity, 
say, the averaged group sunspot number R g (Hoyt et al. 
1994) 1 . This fitting is also shown in Fig. 3 as follows. 
The dashed-dotted line shows the time dependence of the 
latitude-averaged function —6{B 2 } for the same parame- 
ters. The fitting factor 6 is determined by the condition 
6(B 2 ) t=3 . 5 27t d = 0.07. The values ~R g /2l60 are shown in 
Fig. 3 by the filled diamonds. The fitting factor 1/2160 is 
determined by the condition: R g (t = 1988)/2160 = 0.07. 
(Our slightly awkward looking choice of scale for some of 



1 The data for R g were taken from the URL 
ftp.ngdc.noaa.gov/STP/ 
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the quantities plotted in Fig. 3 arises because we choose 
to fit our model in the Northern hemisphere where (x c ) 
and (H c ) are negative.) 

It can be seen from Fig. 3 that the maxima of the 
magnetic energy (B 2 ) nearly coincide with the minima of 
the function \ c - A similar behaviour was observed also 
for the quantity R g (a tracer of cyclic solar activity) and 
H c (the tracer of the current helicity), cf. Bao & Zhang 
(1998). 

6. Discussion and conclusions 

We have presented above a simple model of the nonlinear 
solar cycle with a dynamo saturation mechanism based 
on magnetic helicity conservation. The model is obviously 
oversimplified, however it reproduces some features of real 
stellar cycles as far as they are known from, say, the 
Wilson sample data (Baliunas et al. 1995). In particu- 
lar, we obtain both singly periodic and doubly periodic 
cycles. As expected from qualitative arguments as well 
from observational data (see e.g. Bruevich et al. 2001) the 
transition from singly periodic to doubly periodic cycles is 
associated with a general trend from cycles with smaller 
amplitude to those with larger amplitude. 

Our model of dynamo saturation gives stable oscilla- 
tions for a limited region of parameter space. If the dy- 
namo generation becomes stronger the numerical solution 
runs away. On one hand, we have not included in the 
model all kinds of nonlinear dynamo saturation; e.g. buoy- 
ancy and spot formation obviously lead to some losses of 
toroidal magnetic field and thus contribute to dynamo sat- 
uration. In addition, increasing the dynamo number can 
reduce the radial spatial scale of the toroidal magnetic 
field and effectively enlarge the parameter \x which deter- 
mines the toroidal magnetic field dissipation. If /i is inde- 
pendent of D, we thus artificially overestimate the genera- 
tion effect for larger dynamo numbers. On the other hand, 
experience from dynamo simulations as well as observa- 
tional data suggests that cyclic behaviour is typical for 
moderate dynamo action only and chaotic behaviour oc- 
curs when the dynamo action is stronger. It is more than 
natural to expect that a parameter range with chaotic 
temporal behaviour will also exist in our model. 

We have concentrated our attention here on current 
helicity data, so a deeper comparison of the cycle param- 
eters with observational data is beyond the scope of this 
paper. Also, we should keep in mind that we compare the 
simulated current helicity oc \ c w hh just one part of the 
surface current helicity (6 Z (V x b) z ). 

The properties of current helicity simulated with our 
model have been compared with the available observa- 
tional data. The results of this comparison look quite 
promising in spite of the quite limited extent of the ob- 
servational data, as well as the crude nature of the model. 
We have been able to choose a set of governing parameters 
which give helicity properties comparable with the avail- 
able phenomenology. From another viewpoint, the obser- 
vational data are rich enough to indicate a disagreement 



between the available observations and the predictions of 
the model with other parameter sets. 

The parameters which gives an agreement between 
simulations and observations are quite plausible. However 
we feel that it is too early to insist that this agreement 
is more than a coincidence. Really, we base our compar- 
ison on the 10 year observational data of one scientific 
team. We stress that an extension of the observational 
programme to cover several cycles as well the inclusion 
of data obtained by other scientific teams and from other 
tracers would be very important. In particular it would be 
valuable to include the cross-helicity data into the analy- 
sis. 

Note, that although the available observations cover 
just the period of 10 years, these data extend over parts 
of two different solar cycles, namely the 22nd (1988-1995) 
and the 23rd (1996-1997). The available data suggest that 
the shape of the current helicity distribution is the same 
for both cycles, and further, more recent, current helicity 
studies (Bao et al. 2000; 2002) support this interpretation. 
In particular, the hemispheric rule is obeyed from cycle to 
cycle (e.g. Pevtsov et al. 2001). 

We stress that the scenario described in the present 
paper does not include all possible types of nonlinear 
processes which can occur at the nonlinear stage of the 
dynamo (see, e.g., Brandenburg & Subramanian 2000; 
Brandenburg & Dobler 2001), but rather is restricted to a 
minimal number of processes involved in magnetic helicity 
conservation. In the spirit of the basic ideas about the non- 
linear saturation of solar dynamos, the analysis presented 
here has been restricted to the evolution of a only while 
detailed simulations (e.g. Blackmail & Brandenburg 2002; 
Brandenburg & Sokoloff 2002) also demonstrate a quench- 
ing of the turbulent magnetic diffusivity. A quantitative 
model for a nonlinear quenching of turbulent magnetic 
diffusivity has been recently suggested by Rogachevskii & 
Kleeorin (2001) and used in a galactic dynamo model by 
Kleeorin et. al. (2003). 

We appreciate that alternative interpretations for the 
current helicity observations could be suggested, e.g. he- 
licity could in principle be produced during active region 
formation rather than being connected with the magnetic 
field properties in the region of intensive dynamo action. 
A development of such an alternative explanation to the 
point where it could be confronted with the observational 
data looks highly desirable. 

We conclude with the expression of a guarded but real 
optimism that magnetic helicity observations can result 
in the foreseeable future in a new level of understanding 
in dynamo theory which will base the accepted a-profilcs 
not only on order-of-magnitude arguments and numerical 
simulations, but also on the observational data. It would 
be very important to support this progress in solar dy- 
namo theory by similar progress in other areas of dynamo 
theory, in particular for galactic dynamos. We mention in 
this respect a recent suggestion of Enfilin & Vogt (2003) 
concerning the possibility of observing the magnetic helic- 
ity of galactic magnetic fields. 
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Appendix A: The nonlinear system of dynamo 
equations 

The system of nonlinear equations in nondimensional form 
for an axisymmetric mean magnetic field B — Be^ + 
Vx(j4e0) reads 



dB 

— = D{lA + A s B , 
at 

— = a B + A S A , 
at 



(A.l) 
(A.2) 



dx c , X c _ {2R\ i sin 2 0(VsA) (VsB) _ BAgA 



dt T \ I 
C 
R 



-aB 2 + ^ V • {l 2 K- l x v </> v B 2 e r ) ) + «A X C , (A.3) 



where a = X v 0v + x c <l>m/p(r), A s = sm 2 6V 2 s = A 
l/(r 2 sin 2 6>), and 



By 

V S A 

ClA 



(sin 9 A) , B e 



R a 8 
r sin 9 89 

_j_(e r JL( r A) 
r sm 9 \ or 

1 d(fl,Ar sin 9) 

r 



R a 8 
r dr 



(rA) , 



ee 



1 8 
sin 9 89 



(sm0A) 



8(r, 9) 



Now we consider dynamo action in a thin convective shell, 
average A, B and a over the depth of the shell and con- 
sider these quantities as functions of colatitude 9 only. 
Then we neglect the convection shell curvature and replace 
it by a flat slab. This implies that A s = A = 8 2 /89 2 -fi 2 1 
sm9V s = e e (8/89) + ^e r and ClA = g{8A/89). We 
also redefine C fj,l 2 /A p R as C. This yields Eqs. JSJ, (JSJ 
and fity. 
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